function [USPAI] = USPArecon(bin_file_path, acq_parameter_file_path, save_path, probe_param, modality, recon_method, log_compress_parameter, sound_speed)
% 这是一个提供了各种扩展可能性的US和PA图像重建的函数接口。probe_param和acq_parameter_file分别提供了探头和采集方法的参数
% Boyi RENEW 20231120 for 2.5MHz org US and PA and 10.0MHz US with EPSRI and PA
% bin_file_path            : 原始数据bin文件路径，每个bin文件第一帧为光声信号，其余帧为超声信号，此bin文件需要在前端界面选择
% acq_parameter_file_path  : 原始数据采集参数文件（.csv文件）的路径，前期用不到可忽略，无需在界面选择
% save_path                : 保存图像及图像相关参数mat文件的路径，默认保存在bin_file相同目录下，无需在界面选择
%probe_param               ：探头参数，探头参数再加上acq_parameter
% modality                 : 成像模态，0为超声，1为光声
% recon_method             : 重建算法的方法，默认为DAS，后续可能增加新的重建算法，可在界面选择
% log_compress_parameter   : 对数压缩因子，默认为14，可在界面设定
% sound_speed              : 声速，默认为1500，可在界面设定
% USPAI = USPArecon('E:\Projects\Test\Rawdata_FPGA1_Frame0.bin',1, 'E:\Projects\Test\', 0,'DAS', 16, 1500);

% disp(probe_param.Nelements);
% disp(probe_param.name);
% disp(probe_param.fc);
% disp(probe_param.pitch);
% disp(probe_param.kerf);
% disp(probe_param.width);
% disp(probe_param.bandwidth);
% disp(probe_param.radius);
% disp(probe_param.focus);
% disp(probe_param.height);

[Data_Rf, acqLen] = DataExtractFromBin(bin_file_path, 0);

% (60mm deeep 3200point 1500m/s)    %(45mm deeep 2400point in water 1500m/s)
recon_depth = 2400;
if acqLen < recon_depth
    recon_depth = acqLen;
end

if modality == 0  %NOTE:判US=0/PA=1   
    switch recon_method
        case 'DAS'%NOTE FOR_2.5MHz_HR-DAS
            params.fs = 40e6; %[Hz]
            params.fc = 2.5e6; %[Hz]
            params.bandwidth = 10; % [percent]
            params.pitch = 0.300e-3; % pitch == spacing [mm]
            params.width = 0.254e-3; 
            params.radius = Inf;
            params.c = sound_speed; % [m/s], the reconstruction sound speed - water 1480 m/s
            params.fnumber = []; % determine automatically, for formula see MUST user guide
            % generate reconstruction grid, here set to match imaging area
            xlen = params.pitch *63  + params.width;
            x = (-xlen/2):0.05e-3:(xlen/2);  % spacing = 0.05e-3
            z = 0:0.05e-3:45e-3;% exp：60mm = 60e-3       exp：45mm = 45e-3
            [X,Z] = meshgrid(x,z); % [mm]   
            fs = 40e6;
            fl = 2.1e6;
            fh = 3.0e6;
            % wp 设置截止频率（带通或者带阻） 
            wp=[fl/(fs/2) fh/(fs/2)];
            N=32; 
            % b=fir1(n,Wn,Window ) n滤波阶数，Wn截止频率，0≤Wn≤1，Wn=1对应于采样频率的一半。
            % 当设计带通和带阻滤波器时， Wn=[W1 W2],W1≤ω≤W2。
            % Window — 窗函数。窗函数的长度应等于FIR滤波器系数个数，即阶数 n+1。
            % b返回滤波器系数
            b = fir1(N,wp,blackman(N+1));  %  blackman(n)产生一个长度为n的布拉克曼窗 

            Data_US_Rf = Data_Rf.US(1:recon_depth,:,:); 
            Data_US_Rf=[Data_US_Rf(1:end,:,1)];

            pulser_width = 0.2; %μs @2.5MHz 1-N cycle or 1-element code
            n_cycles     = 2;    
            correct_pulser_point = uint32((pulser_width * n_cycles)*(1e-6)/(1/40e6));%补偿到开始采集点数:时间/采样频率
            Data_US_Rf = [zeros(correct_pulser_point,64,1);Data_US_Rf]; 
        
            Data_US_Rf = filtfilt(b,1,Data_US_Rf); 
            Data_US_Rf = Data_US_Rf-mean(Data_US_Rf,1); 
                    
            IQ_US = rf2iq([medfilt1(Data_US_Rf(1:recon_depth,:,1))],40e6,2.5e6);
            params_US = params;
            bfsig_US = das(IQ_US,X,Z,zeros(1,64),params_US);

%             if ~exist("MUS2_5MHZ","var")
%                 MUS2_5MHZ = load("MUS2_5MHZ.mat");
%                 MUS2_5MHZ = MUS2_5MHZ.MUS;
%             end
% 
%             [IQ_US_m,IQ_US_n]=size(IQ_US);
%             IQ_US_mn=IQ_US_m*IQ_US_n;
%             bfsig_US_temp = MUS2_5MHZ*reshape(IQ_US,[IQ_US_mn,1]);
%            bfsig_US = reshape(bfsig_US_temp,size(X));

            origin_recon_us_image = bmode(abs(bfsig_US),log_compress_parameter);
            if ~exist('USPSF','var')
                load("USPSF.mat", "USPSF");
            end
            WTUS = zeros(size(origin_recon_us_image));
            WTUS(3:end-1,3:end-1) = 1;
            INITUSPSF = USPSF;
            [USJ USP] = deconvblind(origin_recon_us_image, INITUSPSF, 10, uint8(0), WTUS);
            USJ(1:200,:)=0;
            USPAI.image = USJ;
            USPAI.Colormap = colormap("gray");

%             USPAI.Colormap = transpose(colormap("gray"));
%             %USPAI.image = rot90(origin_recon_us_image);
%             %USPAI.image = origin_recon_us_image;
%             USJ(1:200,:)=0;
%             USPAI.image = rot90(USJ);

            set(gca, 'box', 'off',...
                'linewidth',4.5,...
                'TickDir','in',...
                'fontsize',16,'fontname','Times','FontWeight','bold','dataaspectratio',[1,1,1], 'XTick', [1 96 192 288 384],'XTickLabel', [0 5 10 15 20],'YTick', [1 100 200 300 400 500 600 700 800 900],'YTickLabel', [0 5 10 15 20 25 30 35 40 45],'box','on');

            USPAI.dimension = size(USPAI.image);
            USPAI.pixel_spacing = [0.05,0.05];
            USPAI.channels = 1; %获取像素通道数
            USPAI.datatype = class(USPAI.image);
            USPAI.XTick = [0 5 10 15 20];
            USPAI.XTickLabel = '[mm]';
            USPAI.YTick = [0 5 10 15 20 25 30 35 40 45];
            USPAI.YTickLabel = '[mm]';

        case 'RapidDAS'
            params.fs = 40e6; %[Hz]
            params.fc = 10e6; %[Hz]
            params.bandwidth =5; % [percent]
            params.pitch = 0.150e-3; % pitch == spacing [mm] 深超探头
            params.width = 0.127e-3;  
            params.radius = Inf;
            params.c = sound_speed; % [m/s], the reconstruction sound speed - water 1480 m/shttps://itis.swiss/virtual-population/tissue-properties/database/acoustic-properties/speed-of-sound/
            params.fnumber = []; % determine automatically, for formula see MUST user guide
            % generate reconstruction grid, here set to match imaging area
            xlen = params.pitch *63  + params.width;
            x = (-xlen/2):0.025e-3:(xlen/2);     
            z = 0:0.025e-3:22.5e-3;% exp：60mm = 60e-3       exp：22.5mm = 22.5e-3
            [X,Z] = meshgrid(x,z); % [mm]

            fs =40e6;
            fl = 9.5e6;
            fh = 10.5e6;
            % wp 设置截止频率（带通或者带阻） 
            wp=[fl/(fs/2) fh/(fs/2)];
            N=8; 
            % b=fir1(n,Wn,Window ) n滤波阶数，Wn截止频率，0≤Wn≤1，Wn=1对应于采样频率的一半。
            % 当设计带通和带阻滤波器时， Wn=[W1 W2],W1≤ω≤W2。
            % Window — 窗函数。窗函数的长度应等于FIR滤波器系数个数，即阶数 n+1。
            % b返回滤波器系数
            b=fir1(N,wp,blackman(N+1));  %  blackman(n)产生一个长度为n的布拉克曼窗 

        
            Data_US_Rf = Data_Rf.US(1:recon_depth,:,:); 
            Data_US_Rf=[Data_US_Rf(1:end,:,1)];

            pulser_width = 0.05; %μs @10.0MHz 1-N cycle or 1-element code
            n_cycles     = 2;    
            correct_pulser_point = uint32((pulser_width * n_cycles)*(1e-6)/(1/40e6));%补偿到开始采集点数:时间/采样频率
            Data_US_Rf = [zeros(correct_pulser_point,64,1);Data_US_Rf]; 


            Data_US_Rf = medfilt1(Data_US_Rf(1:2400,:,1));
            Data_US_Rf = filtfilt(b,1,Data_US_Rf); 
            Data_US_Rf = Data_US_Rf-mean(Data_US_Rf,1); 
            
            IQ_US = rf2iq([medfilt1(Data_US_Rf(1:recon_depth,:,1))],40e6,10e6);
            params_US = params;
%             bfsig_US = das(IQ_US,X,Z,zeros(1,64),params_US);

            if ~exist('MUS10MHZ','var')
                load("MUS10MHZ.mat","MUS10MHZ");
            end

            [IQ_US_m,IQ_US_n]=size(IQ_US);
            IQ_US_mn=IQ_US_m*IQ_US_n;
            bfsig_US_temp = MUS10MHZ*reshape(IQ_US,[IQ_US_mn,1]);
            bfsig_US = reshape(bfsig_US_temp,size(X));


            origin_recon_us_image = bmode(abs(bfsig_US),log_compress_parameter);
            USPAI.Colormap = transpose(colormap("gray"));
            origin_recon_us_image(1:200,:) = 0;
            % EPSRI processing
            IU = mat2gray(origin_recon_us_image);
            % EPF.parameter
            sigma_s = 6;
            sigma_r = 0.6;
            F_rf_U = EPSRI(IU, sigma_s, sigma_r);
            temp_F_rf_U = ((F_rf_U));
            P_F_rf_U=imadjust((temp_F_rf_U),[0.1 0.95],[0 1],1);
            P_F_rf_U=mat2gray(P_F_rf_U);
            % EPF.processing
            USPAI.image = rot90(im2uint8(P_F_rf_U));

            set(gca, 'box', 'off',...
                'linewidth',4.5,...
                'TickDir','in',...
                'fontsize',16,'fontname','Times','FontWeight','bold','dataaspectratio',[1,1,1], 'XTick', [1 96 192 288 384],'XTickLabel', [0 2.5 5.0 7.5 10.0],'YTick', [1 100 200 300 400 500 600 700 800 900],'YTickLabel', [0 2.5 5.0 7.5 10.0 12.5 15 17.5 20 22.5],'box','on');

            USPAI.dimension = size(USPAI.image);
            USPAI.pixel_spacing = [0.025,0.025];
            USPAI.channels = 1; %获取像素通道数
            USPAI.datatype = class(USPAI.image);
            USPAI.XTick = [0 2.5 5.0 7.5 10.0];
            USPAI.XTickLabel = '[mm]';
            USPAI.YTick = [0 2.5 5.0 7.5 10.0 12.5 15 17.5 20 22.5];
            USPAI.YTickLabel = '[mm]';
    end
elseif modality == 1 %NOTE:判US=0/PA=1   
    switch recon_method
        case 'DAS'
            params.fs = 40e6; %[Hz]
            params.fc = 2.5e6; %[Hz]
            params.bandwidth = 10; % [percent]
            params.pitch = 0.300e-3; % pitch == spacing [mm]
            params.width = 0.254e-3; 
            params.radius = Inf;
            params.c = sound_speed; % [m/s], the reconstruction sound speed - water 1480 m/s
            params.fnumber = []; % determine automatically, for formula see MUST user guide
            % generate reconstruction grid, here set to match imaging area
            xlen = params.pitch *63  + params.width;
            x = (-xlen/2):0.05e-3:(xlen/2);  % spacing = 0.05e-3
            z = 0:0.05e-3:45e-3;% exp：60mm = 60e-3       exp：45mm = 45e-3
            [X,Z] = meshgrid(x,z); % [mm]   
            fs = 40e6;
            fl = 2.1e6;
            fh = 3.0e6;
            % wp 设置截止频率（带通或者带阻） 
            wp=[fl/(fs/2) fh/(fs/2)];
            N=32; 
            % b=fir1(n,Wn,Window ) n滤波阶数，Wn截止频率，0≤Wn≤1，Wn=1对应于采样频率的一半。
            % 当设计带通和带阻滤波器时， Wn=[W1 W2],W1≤ω≤W2。
            % Window — 窗函数。窗函数的长度应等于FIR滤波器系数个数，即阶数 n+1。
            % b返回滤波器系数
            b = fir1(N,wp,blackman(N+1));  %  blackman(n)产生一个长度为n的布拉克曼窗 


            Data_PA_Rf = Data_Rf.PA(1:recon_depth,:,:); 
            Data_PA_Rf=[Data_PA_Rf(1:end,:,1)];
        
            Data_PA_Rf = filtfilt(b,1,Data_PA_Rf); 
            Data_PA_Rf = Data_PA_Rf-mean(Data_PA_Rf,1); 
            
            IQ_PA = rf2iq([medfilt1(Data_PA_Rf(1:recon_depth,:,1))],40e6,2.5e6);
            params_PA = params;
            params_PA.passive = true;            
            bfsig_PA = das(IQ_PA,X,Z,zeros(1,64),params_PA);

%             if ~exist('MPA2_5MHZ','var')
%                 MPA2_5MHZ = load("MPA2_5MHZ.mat");
%                 MPA2_5MHZ = MPA2_5MHZ.MPA;
%             end
% 
%             [IQ_PA_m,IQ_PA_n]=size(IQ_PA);
%             IQ_PA_mn=IQ_PA_m*IQ_PA_n;
%             bfsig_PA_temp = MPA2_5MHZ*reshape(IQ_PA,[IQ_PA_mn,1]);
%             bfsig_PA = reshape(bfsig_PA_temp,size(X));
 
            origin_recon_pa_image = bmode(abs(bfsig_PA),log_compress_parameter);
            if ~exist('PAPSF','var')
                load("PAPSF.mat", "PAPSF");
            end
            WTUS = zeros(size(origin_recon_pa_image));
            WTUS(3:end-1,3:end-1) = 1;
            INITPAPSF = PAPSF;
            [PAJ PAP] = deconvblind(origin_recon_pa_image, INITPAPSF, 10, uint8(0), WTUS);
            USPAI.Colormap = colormap("hot");
            USPAI.image = PAJ;


            USPAI.Colormap = transpose(colormap("hot"));
            %origin_recon_pa_image(1:200,:)=0;
            %USPAI.image = rot90(origin_recon_pa_image);
            %USPAI.image = origin_recon_pa_image;
            %PAJ(1:200,:)=0;
            %USPAI.image  = transpose(PAJ);
    
            set(gca, 'box', 'off',...
                'linewidth',4.5,...
                'TickDir','in',...
                'fontsize',16,'fontname','Times','FontWeight','bold','dataaspectratio',[1,1,1], 'XTick', [1 96 192 288 384],'XTickLabel', [0 5 10 15 20],'YTick', [1 100 200 300 400 500 600 700 800 900],'YTickLabel', [0 5 10 15 20 25 30 35 40 45],'box','on');

            USPAI.dimension = size(USPAI.image);
            USPAI.pixel_spacing = [0.05,0.05];
            USPAI.channels = 1; %获取像素通道数
            USPAI.datatype = class(USPAI.image);
            USPAI.XTick = [0 5 10 15 20];
            USPAI.XTickLabel = '[mm]';
            USPAI.YTick = [0 5 10 15 20 25 30 35 40 45];
            USPAI.YTickLabel = '[mm]';


        case 'RapidDAS'
            params.fs = 40e6; %[Hz]
            params.fc = 10e6; %[Hz]
            params.bandwidth =5; % [percent]
            params.pitch = 0.150e-3; % pitch == spacing [mm] 深超探头
            params.width = 0.127e-3;  
            params.radius = Inf;
            params.c = sound_speed; % [m/s], the reconstruction sound speed - water 1480 m/shttps://itis.swiss/virtual-population/tissue-properties/database/acoustic-properties/speed-of-sound/
            params.fnumber = []; % determine automatically, for formula see MUST user guide
            % generate reconstruction grid, here set to match imaging area
            xlen = params.pitch *63  + params.width;
            x = (-xlen/2):0.025e-3:(xlen/2);     
            z = 0:0.025e-3:22.5e-3;% exp：60mm = 60e-3       exp：22.5mm = 22.5e-3
            [X,Z] = meshgrid(x,z); % [mm]

            fs =40e6;
            fl = 9.5e6;
            fh = 10.5e6;
            % wp 设置截止频率（带通或者带阻） 
            wp=[fl/(fs/2) fh/(fs/2)];
            N=8; 
            % b=fir1(n,Wn,Window ) n滤波阶数，Wn截止频率，0≤Wn≤1，Wn=1对应于采样频率的一半。
            % 当设计带通和带阻滤波器时， Wn=[W1 W2],W1≤ω≤W2。
            % Window — 窗函数。窗函数的长度应等于FIR滤波器系数个数，即阶数 n+1。
            % b返回滤波器系数
            b=fir1(N,wp,blackman(N+1));  %  blackman(n)产生一个长度为n的布拉克曼窗 


            Data_PA_Rf = Data_Rf.PA(1:recon_depth,:,:); 
            Data_PA_Rf=[Data_PA_Rf(1:end,:,1)];


            Data_PA_Rf = medfilt1(Data_PA_Rf(1:2400,:,1));
            Data_PA_Rf = filtfilt(b,1,Data_PA_Rf); 
            Data_PA_Rf = Data_PA_Rf-mean(Data_PA_Rf,1); 
            
            IQ_PA = rf2iq([medfilt1(Data_PA_Rf(1:recon_depth,:,1))],40e6,10e6);
            params_PA = params;
            params_PA.passive = true;
%             bfsig_PA = das(IQ_PA,X,Z,zeros(1,64),params_PA);

            if ~exist('MPA10MHZ','var')
                load("MPA10MHZ.mat","MPA10MHZ");
            end
            [IQ_PA_m,IQ_PA_n]=size(IQ_PA);
            IQ_PA_mn=IQ_PA_m*IQ_PA_n;
            bfsig_PA_temp = MPA10MHZ*reshape(IQ_PA,[IQ_PA_mn,1]);
            bfsig_PA = reshape(bfsig_PA_temp,size(X));


            origin_recon_pa_image = bmode(abs(bfsig_PA),log_compress_parameter);
            USPAI.Colormap = transpose(colormap("hot"));
            origin_recon_pa_image(1:200,:)=0;
            USPAI.image = rot90(origin_recon_pa_image);
            set(gca, 'box', 'off',...
                'linewidth',4.5,...
                'TickDir','in',...
                'fontsize',16,'fontname','Times','FontWeight','bold','dataaspectratio',[1,1,1], 'XTick', [1 96 192 288 384],'XTickLabel', [0 2.5 5.0 7.5 10.0],'YTick', [1 100 200 300 400 500 600 700 800 900],'YTickLabel', [0 2.5 5.0 7.5 10.0 12.5 15 17.5 20 22.5],'box','on');

            USPAI.dimension = size(USPAI.image);
            USPAI.pixel_spacing = [0.025,0.025];
            USPAI.channels = 1; %获取像素通道数
            USPAI.datatype = class(USPAI.image);
            USPAI.XTick = [0 2.5 5.0 7.5 10.0];
            USPAI.XTickLabel = '[mm]';
            USPAI.YTick = [0 2.5 5.0 7.5 10.0 12.5 15 17.5 20 22.5];
            USPAI.YTickLabel = '[mm]';
    end

else
    disp("The reconstruction method is not supported now!");
end

close(gcf);

end


